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I  began  this  research  with  the  goal  in  mind  of  creating  a  computer 
program  that  could  be  used  tactically  to  find  the  best  evasive  maneuver  for 
a  satellite  In  geosynchronous  orbit  under  a  given  attack  through  the  use  of 
continuous  thrust.  I  wanted  this  computer  program  to  be  able  to  track  the 
distance  between  the  attacking  satellite  and  the  satellite  being  attacked 
throughout  the  maneuver.  My  intention  was  to  use  a  very  low  thrust  such 
as  that  available  for  use  in  satellite  attitude  control  and  station  keeping.  I 
was  able  to  make  substantial  progress  through  the  use  of  an  IMSL  routine 
known  as  DVCPR. 

The  biggest  problem  was  in  de-sensitizing  the  subroutine,  DVCPR. 
The  extreme  sensitivity  of  this  subroutine  to  co-state  inputs  proved  to  be 
very  time  consuming.  Hopefully,  my  study  will  save  future  students  some  of 
this  painstaking  drudgery. 

I  owe  a  great  deal  of  thanks  to  my  faculty  advisor,  Lt  Col  Joseph  W. 
Widhalm,  for  his  expertise,  confidence,  and  patience  given  to  me  throughout 
this  project.  However,  my  greatest  thanks  goes  to  my  wife,  Carolyn,  for  all 
of  her  invaluable  encouragement  and  understanding. 
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Abstract 


A  satellite  under  attack  by  another  orbiting  body  relying  on  an 
explosive  device  to  effect  a  kill  has  the  problem  of  avoiding  a  volume  of 
•pace  in  which  its  destruction  is  highly  likely.  To  avoid  this  volume,  the 
attackee  could  use  continuous  low  thrust,  such  as  that  provided  by  the 
electric  propulsion  in  attitude  control  thrusters,  if  its  warning  time  and 
orbital  parameters  were  appropriate.  A  model  is  developed  using  optimal 
control  theory  and  is  solved  numerically  for  the  thrust  direction  using  various 
magnitudes  of  thrust.  The  model  progresses  from  a  one  satellite  solution 
(where  the  threat  is  a  fixed  sphere  in  space)  to  a  two  satellite  solution  in 
which  the  distance  between  the  threat  and  the  target  satellite  is  constantly 
updated  (to  avoid  this  same  threat  sphere).  The  results  are  given  for  several 
values  of  thrust,  several  time- of -flights  (time  of  thrusting  for  the 
maneuvered  satellite),  the  optimized  time  of  flight,  and  for  circular  and  non- 
circular  end  orbits 


OPTIMAL  CONTINUOUS  THRUST  ORBITAL  EVASIVE 
MANEUVERS  FROM  GEOSYNCHRONOUS  ORBIT 


I.  Introduction 

Much  research  has  been  done  on  the  astrodynamics  of  orbital  transfer 
and  rendezvous,  but  very  little  of  the  current  literature  addresses  the 
problem  of  avoiding  an  interception  in  space.  With  more  and  more  military 
capabilities  being  placed  into  satellites  (see  Brandt  (6)),  the  means  of 
intercepting  them  are  improving.  It  stands  to  reason  that  the  problem  of 
avoiding  or  defeating  an  interception  also  needs  attention.  Evasive 
maneuvering  is  one  technique  of  protecting  space  assets. 

Most  of  the  current  literature  on  orbital  evasive  maneuvers  is 
concerned  with  impulsive  thrust  (8,  9,  10,  17,  19,  21,  28).  The  purpose  of 
this  thesis  is  to  examine  the  optimal  orbital  maneuver  that  would  be 
required  by  a  satellite  orbiting  at  geosynchronous  altitude  using  a  constant 
magnitude  continuous  low  thrust  propulsion  system. 

Even  the  simplest  orbital  maneuver  may  use  up  the  satellite’s  limited 
supply  of  propellant,  perhaps  shortening  its  useful  life  or  leaving  it  vulnerable 
to  a  second  attack.  It  becomes  obvious  then,  that  an  optimal  solution  is  to 
minimize  the  time  of  propellant  burn.  This  problem  becomes  synonymous 
with  a  minimum  time  optimization  problem. 

The  problem  is  to  determine  the  optimal  thrust  direction  to  transfer  a 
vehicle  from  some  point  on  an  initial  orbit,  to  some  desired  terminal  orbit  in 
minimum  time,  using  a  constant  magnitude  low  thrust  propulsion  system. 


Tk»  problem  is  approached  using  optimal  control  theory  (See  Bryson  (7)). 
The  three*  direction  becomes  the  control  variable  and,  through  the  equations 
of  motion  a  Two  Point  Boundary  Value  Problem  (TPBVP)  is  formulated. 
For  the  desired  result,  the  problem  is  first  developed  as  a  problem  in  which 
tune  ■  fixed  and  the  problem  is  limited  to  two  dimensions.  An  algorithm  is 
developed  to  move  a  satellite  from  geosynchronous  altitude  to  a  maximum 
mdms  circular  orbit  in  a  specific  time.  There  are  two  constraints  used  on 
the  problem  at  this  point;  first  a  tangency  constraint  which  requires  the 
maneuvering  satellite  to  be  tangent  to  a  threat  sphere  at  ip  and  second  a 
oast  ram  t  which  requires  the  final  orbit  to  be  circular.  Then,  since  some 
fuel  in  this  burn  is  being  used  to  circularize  the  final  orbit,  the  final  circular 
orbit  constraint  is  removed  from  the  algorithm.  This  circular  orbit  at  t, 
oustraint  is  omitted  from  all  remaining  runs.  This  results  in  a  greater 
radius  change  in  the  same  amount  of  time.  The  final  part  of  the  algorithm 
j»*  >rpor»u«  two  satellites  being  monitored  at  the  same  time;  one  of  these 
tmng  the  attacker,  and  the  other  being  our  maneuvered  satellite  (the  target 
>r  attaches)  The  problem  is  then  expanded  to  optimize  time;  thus 
v^ompbshiiig  the  intended  task  of  this  problem. 

The  coordinate  system  selected  for  this  study  differs  from  those 
wrterted  by  Smith  (23)  or  Starr  (25).  Where  they  chose  a  polar  coordinate 
system  the  coordinate  system  chosen  for  this  thesis  is  an  Earth -Centered 
Inertial  Cartesian  System  similar  to  Bowman  (4).  The  main  reason  for  this 
sewrtioo  is  that  an  X-Y  system  is  much  easier  to  visualize  than  a  polar 
system  and  a  system  in  which  it  would  be  easy  to  "see"  what  is  going  on 
•  ith  the  thrust  direction  was  desired  for  this  project.  For  ease  in 
mputiag  the  initial  satellite  position  is  chosen  to  lie  on  the  X  axis,  and 
the  problem  is  kept  two  dimensional.  Simply  stated,  the  initial  problem  as 


ONE  SATELLITE 
COORDINATE  SYSTEM 


stated  in  Bryson  (7:66 — ©8)  is:  given  a  constant -thrust  engine,  Th  —  thrust, 
operating  for  a  given  length  of  time,  tp  we  wish  to  find  the  thrust— direction 
history,  u(Q,  to  transfer  a  satellite  from  geosynchronous  orbit  to  the  largest 
possible  circular  orbit.  (See  figure  1) 

This  problem  is  accomplished  with  good  results  for  a  quarter  orbit 
time -of —flight,  and  for  a  four  hour  time-of-flight.  Problems  encountered 
will  be  discussed  in  the  following  sections. 

The  problem  is  then  expanded  by  removing  the  constraint  of  going  to 

the  largest  possible  circular  orbit,  i.e.  the  orbit  is  now  just  being  changed  to 

the  largest  possible  orbit  away  from  the  threat,  not  to  a  circular  orbit; 

however,  the  tangency  constraint  remains.  This  is  where  this  problem  ceases 
to  be  similar  to  Smith  (23)  or  Starr  (25).  Results  are  obtained  using  the 
same  thrust  magnitudes  and  thrust  durations  (time-of-flight)  as  above. 

Now  comes  the  real  objective;  that  of  actually  tracking  the  threat 
during  the  entire  bum  time,  with  respect  to  the  maneuvered  satellite,  and 
optimizing  the  time— of— flight  for  the  maneuver.  This  means  beginning  again 
with  a  new  set  of  equations  (derived  from  the  initial  equations  of  motion)  for 
two  satellites  within  the  same  problem.  Time  minimization  is  accomplished 
by  parameterizing  time  in  the  equations  of  motion  and  applying  optimization 
techniques  to  come  up  with  the  minimum  time-of-flight  for  the  total 
maneuver  to  increase  satellite  orbit  by  a  specific  distance.  Results  are  also 
obtained  for  the  same  t hi  uats  magnitudes  and  thrust  durations  as  before. 

This  thesis  contains  all  the  mathematical  analysis  in  Section  II,  a 
description  of  the  specific  computer  implementation  of  the  algorithm  and  the 
problems  encountered  in  Section  HI,  a  discussion  of  the  results  in  Section  IV, 
and  conclusion  and  recommendations  in  Section  V. 
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In  Section  II,  the  problem  equations  of  motion  are  listed  for  a  fixed 
time,  maximum  radius  increase  algorithm.  As  stated  above,  one  of  the 
initial  constraints  is  then  removed  and  the  problem  is  reworked  to  get  an 
even  greater  radius  increase.  Therefore,  it  becomes  possible  to  get  a  greater 
threat  avoidance  distance  with  this  second  method.  The  problem  is  then 
reformulated  and  expanded  to  a  two  satellite  problem  in  which  the  distance 
between  the  two  satellites  can  be  monitored  throughout  the  maneuver,  (see 
figure  2)  The  problem  of  keeping  track  of  the  threat  satellite  throughout  the 
bum  time  is  also  discussed.  Section  m  discusses,  in  detail,  the  problems 
encountered  in  all  aspects  and  versions  of  the  algorithm,  and  how  success  is 
finally  achieved  for  all  specified  thrust  magnitudes  and  bum  times.  In 
Section  IV,  the  results  are  discussed  for  various  thrust  magnitudes  and  thrust 
durations  with  the  relationship  between  the  two.  The  results  of  the  time 
minimization  problem  are  also  discussed  in  this  section.  Finally,  the 
conclusions  of  this  thesis  are  discussed  along  with  some  of  the  tactics  of 
orbital  evasion  and  some  suggestions  for  further  work  in  this  area. 


5 


•ttackci 

orbit 


The  equations  of  motion  are  developed  from  two- body  motion  utilising 
the  given  coordinate  systems,  (see  figs.  1  and  2)  In  order  to  keep  things  as 
simple  as  possible,  perturbations  are  neglected  and  only  two-body  orbital 
dynamics  are  considered.  The  basic  equations  come  from  Bate  (3)  and 
Baker  (1).  The  only  forces  acting  on  the  satellite  become  thrust  and  gravity. 
The  angle  between  the  thrust  vector  and  the  X  axis  becomes  the  control. 

Single  Satellite  Problem 

For  the  single  satellite  problem  p  is  the  radius  of  the  satellite 
measured  from  the  center  of  the  earth.  In  the  following  equations  p( 
becomes  x,,  py  becomes  x,,  and  so  on.  The  resulting  state  equations  are. 


*\  -  x,  (1) 

x’a  -  x,  (2) 

x'a  *  -nx,/p*  ♦  Th/(m#-m*t)COS(a)  (3) 

x\  *  -px^/p®  +  Th/(m0-m’t)SIN(a)  (4) 


with  initial  conditions  of: 


x,(0)  *  6.622791181  DU 

(5) 

x,(0)  *  0.0  DU 

(6) 

Xj(0)  =  0.0  DU/TU 

(7) 

x4(0)  =  0.3885792278  DU/TU 

(8) 

The  cost  function  for  this  problem  then  becomes 


J  -  P(t,) 

The  Hamiltonian  of  the  system  is,  therefore 

K  »  XjX,  ♦  X,x4  +  X^-iu^/p*  +  Th/(m#-m’t)COS(n))  +  X^-iu^/p* 

♦  Th/(mc-m’t)SIN(n))  (1C 

Minimising  the  Hamiltonian  with  respect  to  n  results  in: 

K%  «  -X,(Th/(m#-m’t))SIN(«)  +  X4(Th/(m0-m’Q)COS(n)  (11) 

Therefore,  the  resulting  control  is: 

TAN(n)  -  X«/X,  (12] 

The  resulting  co-state  equations  will  now  be  given  with  the  following 
representations: 


*tm  *1  N'S  *•  “  (13) 

COS(n)  -  x7/(l+V/x7,)‘/a  (14) 

SIN(n)  -  xj/d+x,8/^*)1/*  (16) 


giving  the  co— state  equations  as: 


X’,  -  -  SxyX,9/^*  -  SXjXjX^p*)  (16) 

X’,  -  14 Vp*  -  SXjXjVf*  -  SxyXjX,/^)  (17) 

x’T  -  -x,  (18 

x',  -  -x,  (19 


The  remaining  boundary  conditions  come  from  the  two  constraints  on 
the  problem  at  ^  and  from  the  necessary  conditions  as  found  from  the  • 
equation.  The  first  constraint  forces  a  tangency  of  the  maneuvering  satellite 
to  an  imaginary  threat  sphere  at  the  final  time;  whereas,  the  second 
implements  a  circular  orbit  at  if 

-  XjXg  +  XjX4  «  0  (20) 

♦3  -  (P/P«f))l'a  -  y{tf)  -  0  (21) 

where: 

v  -  (x,9  +  x,9)1/9  (22) 

♦  -  Ip  I  +  *,(♦,)  +  va(6a)  (23) 

(Note:  X,(V  «  (i  *  1-4)  )  Therefore  the  necessary  conditions 

become 


MV  *  *s(V  "  Xi/P  +  vixs  ~  Vii4l/2/P#/J 

MV  -  x,(L)  -  Xj/p  +  v,x4  -  VjXJi^/p8/9 


(24) 

(26) 


MM  “  MM  •  *i*i  -  v»/v 
MM  -  X,((,)  -  »,X3  -  VjX^V 


(30) 

(27) 


and  the  necessary  conditions  yield  two  more  boundary  conditions  (at  t,): 

XjX,  ♦  x,x,  -  p  -  x^x,  -  x,x4  -  0  (28) 

XjXj+x^x,  -  p  -  p1/2(xTxs+x#x4)/(vpl/»)  -  0  (29) 

These  last  two  boundary  conditions  come  from  algebraic  manipulation 
of  the  X(tf)  equations  when  one  eliminates  the  »’s. 

Two  Satellite  Problem  (see  figure  2) 

The  two  satellite  problem  equations  are  very  similar  to  the  single 
satellite  problem  except  that  they  utilise  the  vector  between  the  two  satellites 
as  the  position  state  vector  instead  of  a  single  satellite’s  position  vector. 

This  is  done  in  order  to  observe  what  effect  the  dynamics  of  the  attacking 
satellite  has  on  the  maneuvering  problem.  The  following  representations  are 
used: 


e  "  h  ~  Ci. 


e’  -  Ya  -  V. 


*i  *  K  ^  *•  "  P’«> 


The  resulting  state  equations  are: 


(30) 

(31) 


x\  =  x,  (32) 

x’j  -  x4  (33) 

x’,  *=  ~l*n/r2*  +  N’l./r,*  -  Th/(m0-m’l)COS(a)  (34) 

x’4  -  -ju-jy/r a*  +  pj^/r,*  -  Th/(m#-m’l)SIN(a)  (35) 


where  i3  wu  found  through  the  use  of  a  Kepler  orbit  determiner  and 
[,  *u  then  found  by  subtracting  the  state  variables  from  the  appropriate 
components  of  the  ra  vector,  i.e.: 

Ci  ■  (r*  -  Xj)i  +  (r^  -  Xj)i  (30] 

The  following  become  the  final  problem’s  initial  conditions: 


xt(0)  -  -7.652551  DU  (i 

x,(0)  -  -.7744873  DU  (i 

x,(0)  -  1.0241614  DU/TU  (J 

x4(0)  -  -.9642712  DU/TU  (' 


(Note:  These  are  the  result  of  satellite  I  being  in  a  geosynchronous  orbit 
and  satellite  1  being  in  an  eccentric  orbit  with  its  perigee  altitude  '100  nn 
and  its  apogee  altitude  being  geosynchronous.) 

The  cost  function  of  this  system  remains: 

J  -  P(t») 

and  the  Hamiltonian  of  this  system  is,  therefore, 

K  -  XjX,  +  kfn  +y-HJ,ta/r2*  +  — Th/(m0-m*f)COS(«)) 

+  M*"l4ra(r/ra'  +  l“Vri#  -Th/(m,-m,f)SIN(«))  (41) 
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Minimising  this  Hamiltonian  with  respect  to  n  results  in  the  same 


expression  as  before: 

-  X8(Th/(m0-m't))SIN(n)  -  X4(Th/(m0-m*Q)COS(n)  (11) 

Therefore,  the  control  remains: 

TAN(n)  *  X^X,  (12) 

The  co— state  equations  will  remain  given  with  the  following 
representations: 

*s  “  *i  **  “  x,  -  X,  x,«X4  (13) 

COS(n)  -  VUV/x y*)'*9  04) 

SIN(n)  -  x^d+x.Vxr*)1'2  (15) 

giving  the  co-state  equations  as: 

x',  =  rtxj/r,*  “  ^uVr,®  -  axjr^^/r,8]  (42) 

x’*  *  “  3xsri//ri&  -  3x7rltrly/r,#)  (43) 

x’t  =  -xs  (44) 

x’,  *  -x,  (45) 

The  remaining  boundary  conditions  come  from  the  constraint  on  the 
problem  at  t,  and  from  the  necessary  conditions.  Once  again,  the  constraint 
forces  a  tangency  to  the  threat  sphere  at 


where: 


♦  -  P  +  MM  (*« 

the  necessary  conditions  become 


MV  “  MV  “  Vp  ♦  vix« 

(47) 

MV  “  MV  *  Vp  +  \x* 

(48) 

MV  *  MV  “  *ixi 

(49) 

MV  *  MV  * 

(50) 

Sind  the  necessary  conditions  yield  the  remaining  three  boundary 
conditions  (at  (,): 


*!*• 

-  ¥i*° 

(61) 

V 

+  Vs  “  0 

(62) 

Vi 

+  x,Xj  -  p  *  0 

(53) 

Once  again,  these  last  three  boundary  conditions  come  from  algebraic 


manipulation  of  the  k{t^)  equations  when  one  eliminates  v, 


I 


*T* 


W 


To  add  validity  to  the  result*  of  this  thesis,  this  section  begins  with  a 
description  of  how  the  algorithm  actually  works  and  how  it  was  verified. 
Then  the  specific  problem  areas  encountered  in  the  running  of  all  the 
different  algorithm  versions  using  the  different  constraints  are  discussed. 


Algorithm  and  Verification 

The  algorithm  is  written  (all  versions)  in  FORTRAN  77.  It  utilizes 
the  IMSL  routine  DYCPR  to  solve  the  problems  in  question.  It  was  run  on 
the  VAX  11/785  Scientific  Support  Computer  (SSC)  under  the  Unix 
operating  system  at  the  School  of  Engineering  of  the  Air  Force  Institute  of 
Technology  The  algorithm  requires  the  first  derivatives  of  the  equations  of 
motion  as  well  as  the  first  derivatives  of  the  co- state  variable  equations. 

The  first  derivatives  of  the  equations  of  motion  are  multiplied  by  Lagrange 
multipliers  to  form  the  Hamiltonian  of  the  system  The  first  derivatives  of 
the  co  state  variables  (the  Lagrange  multipliers)  come  directly  from  this 
Hamiltonian  equation  as  stated  in  Bryson  (7).  Since  there  are  four  equations 
of  motion  involved  in  the  first  part  of  the  problem  and  four  co-state 
variables  that  come  from  the  Hamiltonian  of  the  system,  the  algorithm 
utilizes  eight  equations  as  input  to  the  IMSL  routine.  These  equations  are 
found  in  the  subroutine  FCNI  This  results  in  an  8  x  8  Jacobian  matrix 
(tx,/^  i*l,8  j*=  1 ,8)  This  Jacobian  matrix  is  found  in  the  subroutine 

FCNJ  The  system  boundary  conditions  are  utilized  in  the  subroutine 


U 


The  basic  discretization  used  is  the  trapezoidal  role  over  a  possibly 
non-uniform  mesh.  This  mesh  is  chosen  adaptively,  in  order  to  make  the 


local  error  of  approximately  the  same  size  everywhere.  Higher  order 
discretizations  are  obtained  by  deferred  corrections.  Global  error  estimates 
are  produced  in  order  to  control  the  computation.  The  resulting  nonlinear 
algebraic  system  is  solved  by  Newton’s  method,  with  step  control,  while  the 
linearized  sparse  system  is  solved  by  a  special  form  of  Gaussian  elimination 
that  preserves  the  sparseness. 

The  algorithm  was  verified  by  inputting  the  required  data  for  running 
a  simple  two  body  problem  with  the  boundary  conditions  split  between  the 
initial  and  final  time.  This  problem  was  run  without  thrust  for  six  hours 
duration  from  geosynchronous  altitudes  and  velocities  and  the  output  resulted 
in  exactly  a  two  body  Keplerian  orbit.  The  final  stage  of  this  thesis 
incorporates  a  Kepler  orbit  determiner  to  find  the  position  of  the  attacking 
satellite  during  the  entire  thrust  duration.  After  verification  of  this  orbit 
determiner,  (through  the  use  of  various  examples  from  Bate  (1:195-210))  it 
was  also  used  to  verify  the  two  body  results  obtained  through  this  algorithm. 

Single  Precision  Version 

Equations  (1)  through  (4)  and  (10)  through  (19)  were  used  as  the 
XPRIME  equations  for  the  subroutine  FCNI.  This  resulted  in  an  8  x  8 
Jacobian  matrix  that  was  used  in  the  subroutine  FCNJ.  Equations  (5) 
through  (8),  (20),  (21),  (28),  and  (29)  were  used  in  the  subroutine  F  JNB  as 
the  boundary  conditions. 

Several  problems  arose  in  obtaining  numerical  results  using  a  single 
precision  version  of  DVCPR.  First,  the  initial  definition  of  the  control  was 
giving  divergent  solutions.  The  computer  was  not  giving  the  correct  values 


because  the  state  equations  needed  COS(u)  and  SIN(u)  inputs.  The  initial 
definitions  were: 

OOS(n)  -  ^/(x,2  +  xf)'/*  (14) 

SIN(u)  -  V(V  +  *7*)l/ *  (15) 

While  mathematically  correct,  this  gave  some  problems  in  the  computer 
algorithm  as  they  caused  the  algorithm  solution  to  diverge  when  first 
attempted.  This  was  rectified  in  two  ways.  First,  i  was  redefined  as 
follows: 


a  *=  ATANS^)  (64) 

This  utilized  the  intrinsic  ATAN2  function  of  FORTRAN  77.  The 
other  way  was  to  painstakingly  take  care  in  choosing  the  input  grid  values  of 
the  co-state  variables.  It  is  interesting  to  note  that  these  two  methods 
resulted  in  exactly  the  same  results;  however,  the  first  solution  was  far  less 
time  consuming. 

The  second  and  biggest  problem  in  using  DVCPR  for  this  problem 
was  that  I  had  to  "initialize  the  grid"  or  give  the  Newton  algorithm  a  good 
starting  point.  Picking  good  approximations  for  the  Initial  state  equations 
was  not  hard,  but  getting  a  good  guess  for  the  co-state  values  proved  very 
time  consuming,  as  DVCPR  proved  to  be  very  sensitive  to  these  Initial 
guesses.  (This  was  the  most  time  consuming  part  of  this  entire  thesis.) 

This  problem  was  finally  rectified  by  utilizing  many  computer  runs  with 
various  co-state  variable  "guesses"  for  the  initial  grid  variables.  This  was 
repeated  until  the  output  began  "looking"  close  to  the  desired  results,  as 


predicted  from  analysis  of  output  from  Smith  (23)  and  Bowman(4).  It  was 
then  possible  to  develop  better  "guesses"  based  on  the  previous  output  and 
sort  of  slowly  "walk"  toward  a  guess  from  which  the  algorithm  would 
converge.  Once  there  was  a  successful  (i.e.  converged)  solution,  it  would  be 
used  as  the  initial  guess  for  subsequent  runs. 

The  third  problem  was  that  of  convergence.  DVCPR  also  proved  to 
be  very  sensitive  to  the  TOL  (tolerance)  parameter.  Small  values  of 
tolerance  in  DVCPR  caused  an  error  message  indicating  that  the  Newton 
algorithm  diverged,  even  though  the  solution  appeared  completely  reasonable. 
Playing  around  with  this  parameter,  i.e.  raising  TOL  from  10~*  to  10 
satisfied  this  convergence  problem  without  affecting  the  actual  output 
tolerances  too  much.  No  explanation  can  be  given  for  this  unusual  behavior, 
as  DVCPR  does  not  document  the  tolerance  algorithm  in  sufficient  detail; 
however,  this  became  one  of  the  primary  reasons  for  developing  the  double 
precision  algorithms. 

Another  problem  was  a  lack  of  convergence  for  any  maneuver  (time- 
of-flight)  that  was  over  approximately  six  and  one-half  hours  in  length. 
Utilization  of  only  one  quarter  (1/4)  of  an  orbit  for  the  runs  avoided  this 
problem.  Justification  for  this  thinking  is  that  in  an  actual  ASAT  type  of 
defense  (in  the  future),  probable  warning  time  will  be  between  four  and  six 
hours  (see  Bowman  (4),  Brandt  (6),  Wagner  (27),  or  Zazworsky  (30))  One 
quarter  of  a  geosynchronous  orbit  is  six  hours,  therefore  six  hours  was  taken 
as  the  maximum  amount  of  time  to  maneuver  in  this  problem. 

The  final  problem  encountered  was  one  of  checking  the  results.  In 
order  to  do  this,  results  were  compared  with  those  obtained  by  Bowman  (4). 
Their  results  utilized  thrust  only  in  the  tangential  direction  and  their 
coordinate  system  was  slightly  different;  however,  the  results  obtained  here 
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compare  favorably  with  their  study.  Result*  also  compare  very  well  with 
those  obtained  by  Smith  (23)  and  Starr  (25). 

Double  Precision  Version 

It  became  obvious  very  quickly  that  single  precision  results  were 
Inconsistent  and  Incomplete,  (see  figures  3-6)  Therefore  a  double  precision 
version  of  the  algorithm  was  developed  for  all  remaining  work.  This  seemed 
to  rectify  most  of  the  above  mentioned  problems  and  orbit  radius  was 
established  to  a  10~*  meter  tolerance.  The  biggest  problem  remained 
however;  "initializing  the  input  grid"  for  the  co-state  values  still  proved  very 
time  consuming.  Once  again,  guessing  until  the  results  were  close  and  then 
using  very  small  variations  of  the  co-state  variables  soon  produced 
convergent  solutions.  The  results  proved  to  be  consistent  with  those 
obtained  by  Bowman  (4)  and  Smith  (23).  The  algorithm  was  initially  run 
with  the  tangency  constraint  and  the  circular  final  orbit  constraint.  (See 
equations  20  and  21)  The  circular  orbit  constraint  was  then  dropped  to 
allow  for  a  larger  (greater  radius)  orbit  in  the  same  amount  of  burn  time. 

(i.e.  more  distance  with  no  more  fuel  than  before)  Results  for  both  circular 
and  non— circular  final  orbits  are  discussed.  Boundary  conditions  are  given  in 
the  subroutine  FCNB  and  are  included  in  the  appendix. 

Two  Satellite  Problem 

Up  until  now,  this  thesis  has  dealt  with  increasing  a  single  satellite’s 
radius  to  the  maximum  attainable  in  a  given  time  subject  to  the  constraints 
of  equations  20  and  21  and  later  dropping  the  constraint  of  equation  21. 
However,  what  this  thesis  was  intended  to  do  is  to  increase  or  maximize  a 
distance  away  from  a  given  threat  satellite  in  an  optimum  time.  This  means 
that  we  must  introduce  this  threat  satellite  into  the  problem,  and  somehow 
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begin  determining  just  how  far  away  from  it  our  own  satellite  is.  At  any 
given  time,  as  we  maneuver  our  satellite  to  avoid  the  threat,  it  would  be 
nice  to  observe  the  effects  on  the  problem  presented  by  the  attacking 
satellite’s  orbital  dynamics.  This  section  deals  with  two  satellites,  a  threat 
or  attacking  satellite  and  a  maneuvering  or  target  satellite. 

First  a  Keplerian  orbit  determining  algorithm  had  to  be  developed. 

This  orbit  determiner  was  derived  from  Bate  (1:195-206)  and  incorporated 
into  the  main  algorithm  in  every  instance  where  the  state  variable  was  to  be 
found.  This  orbit  determiner  found  £2  at  each  time  step,  and  was  then 
found  by  subtracting  the  appropriate  state  variable  from  the  proper 
component  of  (See  equation  36)  This  allows  us  to  treat  as  non— time 
varying  in  the  equations  of  motion. 

Problems  encountered  in  this  version  of  the  problem  at  first  seemed  to 
be  insurmountable.  The  algorithm  was  initialised  with  previous  output,  in 
the  hopes  of  a  quick  convergence.  The  state  variable  initialisation  was 
simply  a  matter  of  subtracting  the  previous  state  variable  from  the  output  of 
the  Keplerian  orbit  determiner  subroutine,  (see  equation  25)  However,  the 
co-state  variable  values  follow  no  logical  pattern.  Once  again,  guess  after 
guess  finally  brought  a  converged  solution  which  could  then  be  used  in 
subsequent  runs. 

This  problem  was  first  run  with  the  threat  satellite  in  the  exact  same 
orbit  as  the  maneuvering  satellite.  In  other  words,  the  threat  and  the  target 
begin  at  the  same  point  in  space  in  the  same  orbit.  The  target  satellite 
then  maneuvers  with  its  constant  low  thrust  system,  thereby  changing  its 
orbit.  This  orbit  changes  with  time,  and  the  radius  is  constantly  increasing; 
or  the  distance  from  the  threat  is  constantly  increasing. 
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These  results  were  then  used  to  initialise  farther  rant  in  which  the 
threat  satellite’s  orbit  is  changed  to  various  attack  runs.  The  final  ran  has 
the  threat  in  a  high  eccentricity  orbit  with  a  perigee  altitude  of  100  nautical 
miles  and  an  apogee  altitude  equal  to  geosynchronous. 

Tim*  M1p1m1**tlon  Version 

The  real  problem  statement  was  to  find  the  optimum  time  to  transfer 
evasively  away  from  a  threat.  This  was  accomplished  by  parameterizing 
time  in  the  original  algorithm  and  using  optimal  control  theory  on  the 
equations  of  motion.  The  cost  function  now  changes  from  that  of  equation  0 
to  become 

tf 

J  -  P«,)  +  /  dt  (55) 

o 

where  we  fix  pi*,)  to  a  specific  magnitude. 

This  results  in  an  additional  state  equation: 

x’,  -  0  (66) 


This  comes  about  by  allowing  orbit  transfer  time  to  vary  between  sero 
and  one,  letting: 


t  *  x,  •  T 


(67) 


Because  p(L)  is  fixed,  the  effective  J  becomes, 


In  order  to  implement  this  extra  state  equation,  each  of  the  previous 
state  equations  must  now  be  multiplied  by  the  new  state  variable  (x*)  and 
equation  (50)  becomes  the  ninth  equation  needed  in  the  subroutine  FCNL 
This  also  brings  an  additional  boundary  condition  into  the  problem,  the 
Hamiltonian  of  this  new  system  (now  different  from  equation  9)  is  now  equal 
to  sero. 

K  -  (XjX,  +  AjX4  +  Xjl-iiXj/p*  +  Th/(m#-m’l)COS(a))  +  XJ-nXj/p1 

♦  Th/(m#-m’i)SIN(a))]  •  x,  +  x,  -  0  (58) 

Equation  58  is  included  as  one  of  the  boundary  conditions  in  the 
subroutine  FCNB.  These  boundary  conditions  can  be  seen  in  the  appendix. 
Since  there  is  an  additional  state  equation,  the  Jacobian  is  now  a  9  x  9 
matrix,  and  is  used  by  the  DVCPR  subroutine  FCNJ.  The  problem  of  grid 
initialization  was  still  present,  but  by  utilizing  previous  data  as  initial  guesses 
for  the  initial  grid,  much  time  was  saved  in  this  parameterized  time  problem 
algorithm. 


IV.  MSVLTS 


This  section  discusses  all  of  the  separate  computer  runs  that  were 
accomplished  for  this  problem.  The  first  part  discusses  the  single  precision 
version  of  the  algorithm.  Double  precision  versions  follow  with  both  end 
time  constraints  discussed  separately.  Then  the  two  satellite  problem  and 
the  different  attacker  orbits  that  were  utilised  is  discussed,  and  the  final 
portion  is  devoted  to  the  time  parameterisation  part  of  this  thesis. 

Saak  Pr&kfrn  vsniga 

The  initial  version  of  the  algorithm  was  written  in  single  precision. 
Four  different  runs  were  accomplished,  .OOlg,  thrust  and  .0005ge  thrust,  each 
utilising  a  four  hour  and  a  six  hour  burn  time.  Each  of  these  runs  were 
subject  to  the  constraints  of  equations  30  and  31.  The  boundary  conditions 
for  these  runs  are  those  found  in  equations  5  through  8,  30,  31,  38,  and  30. 
The  results  are  shown  in  figures  3,  4,  5  and  0  showing  plots  of  m  (Thrust 
direction)  vs  time.  As  stated  previously,  comparison  of  these  results  with 
those  obtained  by  Bowman  (4)  is  very  favorable  and  shows  that  for  thrust 
magnitudes  in  the  range  of  .OOOftg,  up  to  .001  g,,  significant  orbit  changes  are 
possible.  This  algorithm  produced  orbit  radius  changes  greater  than  those 
obtained  by  Bowman  (4)  which  was  expected  since  Bowman  dealt  with 
thrust  only  in  the  tangential  direction.  Thrust  levels  less  than  this  have 
little  or  no  effect  on  orbit  radius.  (NOTE:  g,  is  the  gravitational  constant 
for  the  earth  -lDU/TU,-or-33.1408m/s2) 


THRUST  DIRECTION  (a)  v»  TIME  —  .OOlg,  Thru*  —  4  HR  TOF 

Single  Precision,  Single  Satellite 


U  ra  Tim*  4-Br  (OOlj)  ilo(.  pr«c. 


THRUST  DIRECTION  (u)  v»  TIME  —  .00068.  Thru*  —  4  HR  TOF 

Single  Precision,  Single  Satellite 


Orbit  radius  *u  increased  by  '402  kilometers  io  the  OOO&c,  case  with 
a  6  hour  time -of -flight  (1/4  orbit  period)  la  the  OOlg,  cate,  0  hour 
time -of -flight,  orbit  radios  increased  by  *008  kilometers  This  proved  to  be 
the  largest  amount  achievable  with  the  algorithm  as  written  in  single 
precision.  For  the  0005 g^  case  orbit  radius  increased  only  *355  kilometers 
This  was  the  least  amount  of  the  algorithm.  In  the  .OOlg,,  4  hour  case,  the 
increase  was  '494  kilometers. 

From  the  graphical  results,  one  sees  a  lack  of  consistency  in  the  thrust 
direction  vs  the  thrust  time  The  four  hour  graphs  show  a  rapid  change  in 
thrust  direction;  however,  the  change  from  a  positive  thrust  angle  to  a 
negative  thrust  angle  seems  to  follow  no  set  pattern  In  the  six  hour  graphs 
this  transition  seems  to  be  almost  linear  Although  the  radius  increases 
obtained  with  this  version  seem  to  be  valid  the  thrust  direction  angle  results 
proved  to  be  unacceptable.  This  coupled  with  the  before  mentioned  problem 
with  the  TOL  parameter  prompted  a  double  precision  version  of  the 
algorithm  to  be  incorporated  for  all  additional  runs  and  all  further  results  are 
from  the  double  precision  versions. 

Double  Precision  versions 

The  double  precision  versions  of  the  algorithm  alleviated  the  above 
mentioned  problems  with  thrust  direction  angle  irregularities.  Two  versions 
of  the  algorithm  are  included  because  the  first  has  the  constraint  on  the 
problem  of  the  final  orbit  being  circular  (see  equation  21);  whereas  the 
second  version  does  not  require  a  circular  orbit  at  the  end  of  the  burn 
Both  versions  require  the  maneuvering  satellite  to  be  tangent  to  the  threat 
sphere  at  If.  (see  equation  20)  Boundary  conditions  for  these  runs  are  found 
in  equations  5  through  8,  20,  21,  28,  and  29  for  the  circular  orbit  constraint 


runs  ud  equations  5  through  8,  X),  and  51  through  53  for  the  remaining 
run*  (circular  constraint  removed) 


Circular  Orbit  Constraint  at  t, 

The  first  double  precision  version  of  the  algorithm  was  run  with  the 
same  constraints  as  the  single  precision  version, 

♦j  -  x,x,  +  x^x4  -  0  (20) 

-  (1/pd,))1^  -  v(t,)  -  0  (21) 

The  boundary  conditions  are  as  stated  above.  The  problem  was  run 
with  the  same  thrust  magnitudes  and  durations  as  before.  Accuracies  were 
increased  through  the  TOL  parameter  to  10-4  meters,  (the  final  TOL  value 
was  10~ia)  The  results  can  be  seen  graphically  in  figures  7,  8,  9,  and  10. 

It  was  interesting  that  for  the  same  parameters  as  in  the  tingle 
precision  version,  the  radius  increase  values  were  different.  For  the  OOlg, 
thrust,  six  hour  burn,  the  radius  was  increased  *792  kilometers  in  this 
version  ('184  km  more  than  the  single  precision  version).  However  for  the 
0005c,  thrust,  six  hour  burn,  the  radius  increase  was  only  '392  kilometers 
in  this  version  ('70  km  (ess  than  the  single  precision  version).  The  other 
two  cases  had  similar  results;  the  .OOlg,  thrust,  four  hour  burn,  radius 
increase  was  *472  kilometers  ('12  km  less  that  the  single  precision  version), 
and  the  0005 g,  thrust,  four  hour  burn  had  a  radius  increase  of  only  '240 
kilometers  (115  km  less  that  the  sing.  prec.  version).  While  these  values 
proved  to  be  much  more  accurate  than  the  single  precision  version,  they 
were  not  significant  changes.  The  good  news  about  this  version,  as  can  be 
seen  on  the  graphs,  is  that  the  thrust  direction  now  seemed  to  be  consistent. 


THRUST  DIRECTION  (a)  vs  TIME  —  .OOlg,  Thrust  —  0  HR  TOF 
Double  Precision,  Single  Satellite,  Circular  end  Orbit 


U  n  Tim#  fl-Hr  (.OOlj)  doub.  pr*c. 


Direction  (a) 


THRUST  DIRECTION  (u)  vs  TIME  —  .0005*,  Throat  —  6  HR  TOF 
Double  Precision,  Single  Satellite,  Circular  end  Orbit 


U  r*  Tim*  6 -Hr  (.0005g)  doub.  jrr»c 


figure  9. 
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It  appears  that  the  rapid  thrust  direction  change  that  occurs  in  all  cases  is 
dependent  upon  the  total  time  of  maneuvering.  In  the  six  hour  maneuver 
times,  this  rapid  thrust  direction  change  goes  through  180°;  whereas  in  the 
four  hour  maneuver  times,  it  goes  through  0°.  It  becomes  obvious  that 
thrust  direction  history  becomes  dependent  upon  total  thrust  duration. 

Circular  Orbit  Constraint  at  ^  removed 

It  became  obvious  by  going  through  numerous  output  files,  that  much 
of  the  fuel  being  used  to  make  the  final  orbit  circular  was  in  fact  stealing 
from  orbit  radius  increase.  Since  the  objective  was  to  maximize  p  (orbit 
radius),  this  constraint  (equation  21)  was  removed  in  order  to  get  a  greater 
radius  increase;  therefore,  increase  the  distance  the  target  satellite  was  able 
to  get  away  from  the  attacking  satellite.  In  other  words,  this  version  utilizes 
the  same  amount  of  fuel  to  get  significantly  greater  distances  away  from  the 
threat.  The  boundary  conditions  that  changed  for  this  run  were  due  to 
remc/ing  equation  20,  28,  and  29;  while  including  equations  51  -  53. 

These  results  can  be  observed  graphically  in  figures  11,  12,  13,  and  14. 
From  these  figures,  one  observes  that  the  thrust  direction  change  now  always 
goes  through  0°.  Also,  there  is  no  longer  a  rapid  change  from  positive 
thrust  direction  angles  to  negative  thrust  direction  angles.  The  four  hour 
bums  begin  at  approximately  27°  of  thrust  direction  angle  and  begin  a 
moderately  smooth  transition  (still  highly  non-linear)  to  their  maximum 
negative  values  of  about  -120°.  The  six  hour  burns  both  begin  at 
approximately  45°  and  transition  to  a  maximum  negative  value  of  about 
—90°.  Thrust  direction  angle  is  still  seen  as  highly  dependent  upon  the  total 
bum  time;  the  greater  the  total  bum  time,  the  greater  the  initial  values  of 
thrust  direction  angle. 
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THRUST  DIRECTION  (u)  vs  TIME  —  .001*,  Thrust  —  6  HR 
Double  Precision,  Single  Satellite,  Non-Circular  end  Orbit 


THRUST  DIRECTION  (u)  vs  TIME  —  .0008*,  Thru*  —  6  HR  TOF 
Doable  Precision,  Single  Satellite,  Non -Circular  end  Orbit 


U  *»  Tim*  fl-Rr  (  OOOftf )  noo-drc 


THRUST  DIRECTION  (u)  vs  TIME  —  .OOlg.  Thrust  —  4  H] 
Doable  Precision,  Single  Satellite,  Non-Circular  end  Orbit 


The  significant  part  of  this  version  is  the  distance  increase.  Orbit 
radios  was  increased  '1707  kilometers  in  the  .OOlg,  six  hoar  barn.  This  is 
915  kilometers  more  than  the  previous  version!  The  .0005go  six  hoar  barn 
increased  radios  by  '862  kilometers,  almost  500  km  more  than  that 
previously  obtained.  The  four  hoar  barns  were  equally  impressive.  The 
.OOlg,  four  hour  barn  increased  radios  by  *665  kilometers;  while  the  .0005^ 
four  hour  burn  increased  the  radius  by  '335  kilometers. 

Two  Satellite  Algorithm  Results 

The  results  from  the  two  satellite  portion  of  this  algorithm  are  shown 
in  figures  15,  and  16.  Once  again,  the  only  constraint  on  the  two  satellite 
problem  is  the  tangency  requirement  of  equation  20.  Figure  15  shows  both 
satellites  beginning  at  the  same  point  in  space.  The  terminal  boundary 
conditions  are  those  of  equations  20  and  51  -  53.  One  quarter  of  a 
geosynchronous  orbit  is  shown  as  that  of  the  attacking  satellite.  The  target 
satellite  is  the  one  that  has  the  ever  increasing  radius.  It  is  interesting  to 
note  that  not  only  is  the  radius  increasing,  but  the  period  of  its  orbit  is  also 
changing  ever  so  slightly.  This  ever  so  slight  period  increase  has  a 
tremendous  effect  on  this  problem. 

The  threat  satellite's  orbit  was  then  gradually  changed  and  the 
algorithm  continually  updated  with  previous  output  to  keep  a  good 
"initialization  grid”  until  the  attacking  satellite  was  in  a  very  high 
eccentricity  orbit  with  a  radius  of  perigee  of  approximately  100  nautical  miles 
and  an  apogee  altitude  equal  to  geosynchronous.  The  initial  conditions 
eventually  become  those  of  equations  37  -  40.  The  algorithm  was  then 
repeated  for  both  of  the  previous  thrust  magnitudes  with  both  of  the  thrust 
burn  times  with  the  same  terminal  conditions  as  before. 
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Two  Satellite*  Target  initially  in  Geosynchronous 
Attacker  -  High  Eccentricity 


Two  SoUUit*  Owtroll  Picture 


The  result*  from  this  were  nothing  short  of  great!  The  biggest 
increase  in  distance  came  not  only  from  the  increase  in  the  maneuvered 
satellite’s  radius,  as  was  first  assumed,  bat  also  from  the  fact  that  as  the 
radius  was  increasing,  the  satellite  was  also  slowing  down  in  it*  orbit  just 
enough  to  cause  it  to  arrive  at  the  intercept  point  area  after  the  threat 
satellite  had  already  left  this  area.  The  combined  effect  of  increased  period 
and  increased  radius  greatly  increased  the  total  miss  distance  (the  distance 
between  the  attacker  satellite  and  the  target  satellite)  For  the  OOlg,  thrust, 
six  hoar  barn  time,  total  distance  between  the  two  satellites  at  intercept 
time  is  now  *2800  kilometers  The  0006g,  thrust,  six  hoar  barn  tune, 
resulted  in  an  *1413  kilometers  distance  between  the  satellites  The  OOlg, 
thrust,  four  hour  burn  time,  gave  an  approximate  1000  kilometers  and  the 
.OOOfig,  thrust,  four  hour  burn  time,  increased  the  distance  jost  over  500 
kilometers. 

It  is  important  to  note  that  the  thrust  direction  histones  for  the  two 
satellite  runs  matched,  exactly,  the  runs  made  previously  with  the  fixed 
threat  sphere  and  a  single  maneuvering  satellite.  Therefore,  figures  11  -  H 
are  also  valid  representations  of  the  two  satellite  versions  of  this  algorithm 

An  example  of  the  final  two  satellite  transfer  can  be  seen  in  figure  16 
This  figure  shows  that  for  a  six  hour  burn  time,  the  attacker  is  just  before 
perigee  in  its  orbit  at  the  beginning  of  the  barn  The  two  satellites  continue 
in  their  respective  orbits,  with  the  target  satellite  continually  changing  it* 
orbit.  At  the  intercept  time  (six  hoars  in  this  depiction),  when  the  threat 
satellite  is  programmed  to  explode  (or  whatever),  the  target  satellite  is 
approximately  2800  kilometers  further  oat  and  behind  the  threat,  hopefully 
and  apparently  out  of  danger 
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Now  that  the  important  part  of  the  problem  has  been  worked,  it 
leaves  only  the  problem  of  doing  thin  problem  in  the  optimum  time.  Section 
HI  tells  how  this  algorithm  was  accomplished.  A  fixed  orbit  radios  increase 
was  selected,  and  the  minimum  time  to  get  to  that  orbit  radios  was  found. 
This  problem  was  accomplished  with  both  the  .001g,  thrust  magnitude  and 
the  0006gc  thrust  magnitude.  The  problem  was  first  programmed  for  a 
single  satellite  and  then  attempted  for  the  two  satellite  problem.  A  working 
model  of  the  minimum  time  algorithm  for  the  two  satellite  problem  could 
not  be  achieved.  Results  achieved  are  for  a  single  satellite  with  a  fixed 
threat  sphere.  The  tangency  constraint  remains  the  only  constraint  on  this 
problem  and  the  boundary  conditions  are  included  in  the  appendix  for  this 
algorithm.  Figure  17  shows  the  thrust  direction  angle  vs  the  parameterized 
time  for  the  001  g,  thrust  magnitude  case  in  which  target  radius  was 
increased  oy  500  kilometers.  Five  hundred  kilometers  was  used  by  Burk  (8), 
but  is  not  known  to  be  sufficient  for  all  possible  attacks  in  space.  The 
minimum  time  for  a  500  km  threat  sphere  was  found  to  be  12.0085  TU  or  2 
hours,  41  minutes,  28.0  seconds.  Target  thrust  direction  angle  begins  at 
about  21*  and  transitions  fairly  smoothly  toward  a  maximum  negative  value 
of  about  -128*. 

Figure  18  shows  the  thrust  direction  angle  vs  the  parameterized  time 
for  the  0005 g,  thrust  magnitude  case  in  which  target  radius  was  also 
increased  by  500  kilometers.  (500  km  threat  sphere)  The  minimum  time  in 
this  case  was  found  to  be  16.7495  TU  or  3  hours,  45  minutes,  13.7  seconds. 
Here  the  target  thrust  direction  angle  begins  just  above  27*  and  proceeds  to 
a  maximum  negative  value  of  -107*. 


A  run  was  also  made  with  .001  go  thrust  magnitude  in  which  target 
radius  was  increased  by  1000  km.  The  minimum  time  for  this  case  was 
16.7532  TU  or  3  hours,  45  minutes,  16.7  seconds.  In  other  words,  with 
twice  the  thrust  magnitude,  it  was  possible  to  increase  the  radius  twice  as 
far  in  almost  the  same  time.  This  proved  to  be  merely  coincidental 
however,  as  other  runs  did  not  prove  this  linear  relationship  out.  For  a 
1500  km  orbit  radius  increase  with  .001  go  thrust  magnitude,  the  time  was 
25.3598  TU  or  5  hours,  41  minutes,  0.6  seconds. 

The  minimisation  of  the  two  satellite  problem  was  attempted  to 
incorporate  all  aspects  of  the  original  problem  statement;  that  being  to  get 
the  target  satellite  as  far  away  from  the  attacking  satellite  in  the  minimum 
amount  of  time.  The  attacking  satellite’s  orbital  dynamics  coupled  with  the 
complex  dynamics  of  the  minisation  problem  proved  to  be  too  much  for  the 
DVCPR  routine.  Figure  19,  however,  shows  how  this  procedure  would  most 
likely  appear.  In  this  figure,  The  attacking  satellite  is  actually  just  above 
the  X  axis,  on  its  way  to  its  apogee  altitude  which  will  be  coincident  with 
the  target  satellite  at  geosynchronous  altitude.  The  target  satellite;  however, 
begins  its  thrust  maneuver,  and  a  /ew  short  hours  later  is  presumedly  safely 
away  from  the  threat  instead  of  coincident  with  it. 


Minimum  Time  Problem 
.001  Thrust  —  600  km  radius  increase 


U  t>  Time  (.OOtf)  Mis  Tim* 


2.0E-01  4.0E-01  S.0E-01  8.0E-01 


Tim*  *  12.0080 


figure  17. 
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flfure  II. 


Y  component  (DU) 


Proposed  Minimum  Time  Problem 
.OOlg^  Thrust  —  1000  km  threat  distance 


Two  Satellite  Overall  Picture 


figure  19. 
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V,  CfracluglQM.  TiKticg,  and  Suggestions 


Conclusions 

The  problem  of  optimal  time  transfer  to  some  desired  terminal  radius 
increase  using  a  constant  magnitude,  low  thrust  propulsion  system  has  been 
investigated. 

Numerical  solutions  to  this  highly  non-linear  TPBVP  were  attainable 
through  the  IMSL  routine  DVCPR.  Although  this  routine  must  be  initialized 
(rather  accurately)  in  order  to  find  a  converged  solution,  it  seems  that  it 
would  be  applicable  to  a  wide  class  of  orbital  transfer  problems. 

For  geosynchronous  orbits,  it  is  now  determinable  that  the  radius 
increase  has  a  non-linear  relationship  that  is  proportional  to  the  burn  time. 

It  is  obvious  to  say  that  the  longer  the  burn  time,  the  greater  the  orbit 
radius  will  be;  however,  it  becomes  significant  that  this  happens  very  rapidly 
between  the  four  and  six  hour  points.  Over  one  thousand  kilometers  is 
added  to  the  orbit  radius  between  these  times  in  the  .001  g^  thrust  magnitude 
case.  In  the  .0005^  thrust  magnitude  case,  over  five  hundred  kilometers  is 
added  to  the  radius  between  the  four  and  six  hour  points. 

Initial  thrust  direction  angle  also  is  greatly  dependent  upon  the  total 
time  in  question.  The  greater  the  total  burn  time,  the  greater  the  starting 
value  of  the  thrust  direction.  For  the  minimum  time  cases,  this  does  not 
appear  to  be  true  however.  It  appears  that  the  smaller  the  thrust,  the 
greater  the  beginning  value  for  thrust  direction. 

The  key  conclusion  however  is  that  a  constant,  low  thrust  propulsion 
system  does  seem  to  be  a  valid  alternative  to  the  large  impulsive  propulsion 
system  studies  in  the  current  literature  in  accomplishing  evasive  orbital 
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maneuver*.  In  view  of  the  accuracy  of  theae  results,  this  approach  forma  the 
baaia  of  an  effective  numerical  technique  for  completely  general  optimal  low 
throat  orbital  transfer. 

Tactk* 

Many  practical  aapecta  of  the  attack  avoidance  problem  have  not  been 
dealt  with  in  thia  theaia.  In  a  real  attack  on  a  satellite,  the  threat  would 
have  to  be  detected  and  then  tracked  for  a  while  to  determine  ita  orbit. 

From  thia  information  the  defender  could  deduce  the  intended  target  and  the 
time  of  intercept.  The  defender  would  have  to  come  to  aome  conclusion 
regarding  the  lethal  radius  of  the  threat,  either  from  intelligence  information 
or  from  "lucky"  guessing.  He  might  decide  to  defeat  the  threat  by  evasive 
maneuvering  (as  this  thesis  discusses),  by  attacking  the  attacker,  by  attacking 
the  enemy’s  control  system  or  communications,  by  using  chaff  or  decoys,  or 
by  relying  on  his  satellite’s  hardening.  He  could  also  decide  to  do  nothing  at 
all  and  Just  bear  the  loss  of  a  satellite. 

As  was  stated  earlier,  the  earlier  the  maneuver  begins,  the  more 
distance  away  from  the  attacker  the  satellite  can  achieve.  The  attacker  may 
detect  this  maneuver,  and  re-establish  the  intercept  by  a  conn  term  aaeuver 
on  its  part.  This  would  then  have  to  be  detected  and  another  evasive 
maneuver  performed  by  the  target.  Depending  upon  available  fuel  onboard, 
this  could  be  a  losing  effort.  This  all  becomes  a  type  of  cat  and  mouse 
game,  with  each  player  attempting  to  out  guess  his  opponent.  Under  these 
circumstances,  the  orbital  evasion  problem  becomes  the  classical  military 
game  of  trying  to  outfox  the  enemy. 


There  ere  many  other  areas  of  the  orbital  evasion  problem  that  could 
be  addressed,  but  four  more  areas  of  the  problem  as  addressed  in  this  thesis 
will  be  mentioned.  They  are  as  follows:  the  performance  of  the  algorithm, 
the  effects  of  uncertainty,  the  accuracy  of  the  model,  and  the  development  of 
a  three  dimensional  model. 

This  algorithm  proved  to  give  good  results,  but  only  after  much 
frustration  and  time  was  spent  in  finding  good  initial  values  in  which  to 
begin  the  problem.  A  method  should  be  found  to  make  the  algorithm 
reliable  without  having  to  initialise  the  grid  for  the  Newton  algorithm.  One 
alternative  approach  to  this  may  be  to  have  a  library  of  example  data,  for 
various  types  of  orbits  and  maneuvers.  It  was  found,  that  once  a  fairly  good 
set  of  data  was  available  ss  a  starting  guess,  the  algorithm  achieved  superb 
results. 

The  algorithm  in  this  thesis  does  not  account  for  uncertainties  in  the 
input  data.  How  is  it  determined  that  a  specific  satellite  is  under  attack? 
With  the  increasing  density  of  satellites  in  geosynchronous  orbits,  it  will 
remain  very  hard  to  establish  which  satellite  (or  satellites)  needs  to  be 
maneuvered.  Once  this  is  determined,  can  the  required  thrust  magnitude  or 
mors  importantly  the  required  thrust  direction  angle  be  programmed  into  the 
satellite  accurately?  Since  the  thrust  direction  angle  is  constantly  changing, 
some  means  of  communicating  with  the  satellite  must  be  maintained 
throughout  the  maneuver  or  pre-programmed  into  the  satellite  guidance 
system. 

This  algorithm  does  not  take  into  account  any  perturbation  effects  that 
effect  all  orbiting  bodies,  nor  does  it  account  for  a  maneuvering  threat  or  a 
sou -spherical  threat  volume.  Perturbation  algorithms  are  available  but  some 


are  classified  and  others  can  not  be  implemented  with  available  resources.  It 
has,  however,  been  observed  in  the  open  literature  that  a  maneuvering  threat 
is  currently  beyond  the  scope  of  what  our  current  enemies  capabilities  are. 
Also,  any  other  shaped  threat  volume,  would  have  to  be  addressed  in  a 
separate  algorithm. 

The  last  topic  for  further  work  is  obviously  a  three  dimensional  model. 
Thk  thesis  dealt  only  with  geosynchronous  orbits  and  co-planer  attacking 
satellites.  In  order  to  expand  this  algorithm  to  attacking  satellites  in  inclined 
orbits,  a  three  dimensional  model  must  be  developed.  This  model  may  also 
address  the  problem  of  returning  the  satellite  to  its  mission  orbit  in  an 
optimal  or  minimum  fuel  type  of  burn. 


Appendix 


Compqter  Ifrtinei  of  DVCPR  gqbrpotinee  FCNI  FCNJ,  nnd  FCNB 


Smffr  SteUit t  Vtrtio* 


SUBROUTINE  FCNI(N,T,X,XP) 

DOUBLE  PRECISION  T.XM.XPtNLMU.RO.B.D.TH.MO.MD.U 
INTEGER  N 

COMMON  /A/  MU  /B/  TH.MO.MD 
RO-DSQRT(X(i)**2+X(2)"2) 

IF  (RO  .EQ.  0)  RO— 1.0 
B-TH/(M0-MD*T) 

U«DATAN2(X(8),X(7)) 

D»X(7)*X(1)  +  X(8)*X(2) 

XP(1)-X(3) 

XP(2)«X(4) 

XP(3)--MU*X(l)/RO**3  +  B*DCOS(U) 

XP(4)— MU*X(2)/RO”3  +  B«DSIN(U) 
XP(5)-MU*(X(7)/RO**3  -  3.0d+00*X(l)*D/RO**5) 
XP(6)«MU*(X(8)/RO**3  -  3.0d+00*X(2)*D/RO”5) 

XP(7)—  X(5) 

XP(8)—  X(«) 

RETURN 

END 


SUBROUTINE  FCNJ(N,T,X,PD) 

DOUBLE  PRECISION  T.XW.PDIN.NJ.MU.RO.B.C.D.TH.MO, 
t  MD,E,F,G,H,U 
INTEGER  N 

COMMON  /A/  MU  /B/  TH,M0,MD 
RO-D5QRT(X(l)**2+X(2)”2) 

IF  (RO  EQ.  0)  RO-l.O 
B— TH/(M0-  MD*T) 

U-DATAN2(X(8),X(7)) 

F-X(8)/X(7) 

C-(1.0d+00  +  F**2) 

D-X(7)*X(1)  +  X(8)*X(2) 

E-X(8)*X(1)  +  X(7)*X(2) 

G-3.0d+00*X(7)*X(l)  +  X(8)*X(2) 

H-3.0d+00*X(2)*X(8)  -f  X(1)*X(7) 


PD(1,1)«0.0 

PD(l,2)-0.0 

PD(l,3)-1.0d+00 

PD(1,4)«0.0 

PD(l,5)-0.0 

PD(1,«)«0.0 

PD<1,7)«0.0 

PD(l,8)-0.0 

PD{2,1)-00 

PD(2,2)»0.0 

PD(2,3)-0.0 

PD(2,4)*»1.0d+00 

PD(2,5)»0.0 

PD(2,0)«0.0 

PD(2,7)-0.0 

PD(2,8)«0.0 

PD(3, 1)»MU*(-  1.0d+00/RO**3  +  3.0d+00*X(l)**2/RO**5) 

PD(3l2)»3.0d4-00*MU*X(l)*X(2)/RO**5 

PEH3,3)-0.0 

PD(3,4)»0.0 

PD(3,5)-0.0 

PD(3,fl)-0.0 

PD(3,7)«B*DSIN(U)*X(8)/(C*X(7)**2) 

PEK3.8)— B*DSIN(U)/(C*X(7)) 

PD(4,1)-PD(3,2) 

PD(4,2)»MU*(- 1 .0d+00/RO**3  4-  3  0d4-00*X(2)**2/RO**5) 

PD(4,3)-0.0 

PD(4,4)-0.0 

PD(4,5)-0.0 

PD(4,©)»0.0 

PD(4,7)— B*DCOS(U)*X(8)/(C*X(7)**2) 
PD(4,8)«B*DCOS(U)/(C*X(7)) 

PD(5fl)-MU*(-3.0d+00*G/RO**5  +  15.0d4-00*X(l)**2*D 
t  /R0**7) 

PD(5,2)-MU*(-3.0d+00*E/RO**5  4-  15  Od4-00*X(1)*X(2)*D 
I  /RO**7) 

PD(5,3)-0.0 
PEKS^J-O.O 
PD(5,5)« 0.0 
PD<5,6)-0.0 

PD(5,7)"MU*(1.0d4-00/RO**3  -  3.0d4-00*X(l)**2/RO**5) 
PD(5,8)— 3.0d4-00*MU*X(l)*X(2)/RO**5 


PD(«,1)=PD(5,2) 

PD(6,2)»MU*(-3.0d+00*H/RO**5  +  15.0d+00*X(2)**2*D 
t  /R0**7) 

PD(8,3)=0.0 

PD(6,4)=0.0 

PD(6,5)=0.0 

PD(6,6)=0.0 

PD(6,7)=PD(5,8) 

PD(0,8)=MU*(1.0d+00/RO**3  -  3.0d+00*X(2)**2/RO**5) 

PD(7,1)=0.0 

PD(7,2)=0.0 

PD(7,3)=0.0 

PD(7,4)=0.0 

PD(7,5)— 1.0d+00 

PD(7,6)=0.0 

PD(7,7)=*0.0 

PD(7,8)=0.0 

PD(8,1)*00 

PD(8,2)=0.0 

PD(8,3)*0.0 

PD(8,4)=0.0 

PD(8,5)*0.0 

PD(8,e)— 1.0d+00 

PD(8,7)~0.0 

PD(8,8)=0.0 

RETURN 

END 


SUBROUTINE  FCNB(N,XA,XB,F) 

DOUBLE  PRECISION  XAINLXBINLFINJ.KS.KM.RB 
INTEGER  N 

COMMON  /CONV/  KM.KS 
RB->DSQRT(XB{1)**2+XB(2)**2) 

F(1)-XA(1)  -  KM 
F(2)-XA(2) 

F(3)-XA(3) 

F(4)-XA(4)  -  KS 
F(5)«XB(1)*XB(3)+XB(2)*XB(4) 
F(«)-XB(7)*XB(3)+XB(8)*XB(4) 
F(7)-XB(5)*XB(1)+XB(«)*XB(2)  -  RB 
F(8)-XB(8)-XB(7)*XB(2)/XB(1) 

RETURN 

END 


Two  Satellite  Version 


SUBROUTINE  FCNI(N,T,X.XP) 

DOUBLE  PRECISION  T.XINLXPINLMU.B.D.TH.MO.MD.RIH), 
t  R2(4),R2i,R32,R23)R24,V2i(V2a1V23,V24,R20(4)lV20(4), 
t  C,E,F 
INTEGER  N 

COMMON  /A/  MU  /B/  TH,M0,MD 

COMMON  /RS/  R21)R22)R23,R24,V21,V22,V23,V24 

R20(l)-R21 

R20(2)«RM 

R20<3)-R23 

R20(4)-R24 

V30(l)-V31 

V30(2)-V22 

V20<3)-V23 

V30(4)«V24 

CALL  EKPLER(T,MU,R20,V20,R2,V2) 

Rl(l)-R2(l)-X(l) 

Rl(2)-R2(2)-X(2) 

Rl(3)«0.0 

B—TH/(MO-  MD*T) 

C-X(7)*R1(1)  +  X(8)*R1(2) 
D«DSQRT(Rl(l)**2+Rl(2)**2+Rl(3)**2) 
E-DSQRT(R2(1)**2+R2(2)**2+R2I3)**2) 
F-DSQRT(X(7)**2+X(8)**2) 

XP(1)-X(3) 

XP(2)«X(4) 

XP(3)—  MU*R2(1)/E**3  +  MU*R1(1)/D**3  -  B*X(7)/F 
XP(4)— MU*R2(2)/E**3  +  MU*RI(2)/D**3  -  B*X(8)/F 
XP(5)-MU*(X(7)/D**3  -  3.0d+00*Rl(l)*C/D**5) 
XP(«)-MU*(X(8)/D**3  -  3.0d+00*Rl(2)*C/D**5) 

XP(7)—  X(5) 

XP(8)--X(6) 

RETURN 

END 


SUB10UTIN1  rC!*I(N,T,X,PD) 

DOUBLE  PRECISION  T,X(N),PD(N,N),MU,B,C)D,TH,MO,MD,&, 
8  F1G,R3(4),R21,R32,R33>R24,V3l,V23,V231V24.R30(4) 

8  VJ0<4),R1(4) 

INTEGER  N 

COMMON  /A/  MU  /B/  TH.MO.MD 

COMMON  /RS/  R21,R22,R33  R34,V3I  V22.V23.V34 

R20<l)-R21 

R20<2)-R22 

R20(3)-R33 

R20<4)-R24 

V30(l)-V31 

V20(3)-V23 

V20(3)-V23 
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v  - “  A  V-  t-N.V. 


V20(4)«V24 

CALL  KSPL£B(T,MU.R30.V30.R2.V2) 

Rl(l)»R3(l)-X(l) 

Rl(2)-R2(2)-X(2) 

Rl(3  )-00 

B-TH/(MO- MD*T) 

C»X(7)*R1{1)  +  X(8)*Rl(2) 
D-DSQRT(R1(1)**2+Ri(2)**2+Rl(3)**2) 
E-DSQRT(R2(l)**2+R2(2r*2>R2(3)**2) 
F«DSQRT(X(7)~2  +  X(8)**2) 

G-X(7)*R1(2>  +  X(8)*Rl(l) 

PD(1.1)«0  0 
PD<1,2)«00 
PD<1,3)-1  0d+-00 
PD<1.4)-00 
PD<l,5)-00 
PD(l.6)-00 
PD<1,7)~00 
PDO,8)«0  0 

PCH2.1)— 0  0 
PD<2,2)-0  0 
PD<2.3)-0  0 
PD<2,4)-1  Od^OO 
PD(2>3)-0  0 
PCH2,8)-0  0 
PD<2.7)-0  0 
PD<2,8)-0  0 

PD<3,l)-Ml’*(  lOd-OO/D—3  ♦  3  Od  *00*Rllll“  2/ L>**31 

PD<3,2}«3  Od  *00*MVRl|  l  )*RU2)/t>**5 

PD<3,3)-0  0 

PD<3,4)-0  0 

PD(3.5)’-0  0 

PD(3,6)—0  0 

PEH3>7)»  B/F  *  B*\r>**2  F*M 
PD(3.8>»B-\|7r\l§)/F**3 


PD<4  U-PIX3  i) 

PCK4J)»MI'*(  I 
PCK4.3)-0  0 
PDK4)-0  0 
PD<O)-0  0 
PDt4,e)»0  0 
PCH4.7)«»  PtXi.il 

pn<4  8>-  BF  * 


•3d  3  •  3  (3d -00*KM2)#*J  l •  S i 


B*\(8C*J  f  < 


M 


>  >  v  %  v 


PD(5.1)*MlT*((6  Qd-*-00*X(7)*Rl(l)  +  3.0d-t-00*C)/D*’*5 
15  Od  +  00*R1(1)**2*C/D**7) 

PD(5,2)-Ml'*(3  0d-00*G/D**5  +  15  Od-hOO*Rl(l)*Rl(2)*C 
/D**7) 

PD(5,3)=0  0 
PD(5,4)-0  0 
PD<5,5)=0  ... 

PD(5,«)=0  0 

PD(3.7)«MU*(1  Od+OO/D—S  -  3.0d+00*Rl(l)**2/D*-5) 
PD<5.8)=  3  Od+0O*MU*Rl(l)*Rl(2)/D**5 

PD(0,1)*PD(5,2) 

PP<6.2)=MU*((6  0d-+00*X(8)*Rl(2)  +  3.0d+00*C)/D**5 
15  0d  +  00*Rl(2)**2*C/D**7) 

PD(6,3)=0  0 
PD(«,4)-0  0 
PD(6.5)»0  0 
PD(0.0)=00 
PD(C,7)*PD(5,8) 

PD(6  8)-Ml7*(l  Od+O0/D**3  -  3  0d+00*Rl(2)**2/D**5) 

PD(7,1)*0  0 
PD(7,2)~0  0 
PD(7,3)-0  0 
PD<7,4)~0  0 
PD(7,5)-  1  Od  +00 
PD(7,0)-0  0 
PD(7.7)-0  0 
PD(7.8)=0  0 


PD(8,l)-0  0 
PD(8.2)=K)  0 
PD<8, 3)^0  0 
PD<8,4)~0  0 
PD<8.5)»0  0 
PD(8,0)-  1  Od  -00 
PD<8.7)-0  0 
PD(8.8)-0  0 


RETVRN 


SUBROUTINE  rCNB(N,XA,XB,F) 

DOUBLE  PRECISION  XA(N),XB(N),F(N),KS,KM,RB 
INTEGER  N 

COMMON  /CONV/  KM.KS 
RB=DSQRT(XB(1)**2+XB(2)**2) 

F(1)=XA(1)  +  7.552551d+00 
F(2)»XA(2)  +  .7744873d+00 
F(3)=*XA(3)  -  1.0241014d+00 
F(4)*XA(4)  +  ,9G42712d+00 
F(5)=XB(1)*XB(3)+XB(2)*XB(4) 

F(6)=XB(8)  -XB(7)  *XB(2)  /XB(  1 ) 

F(7)=XB(3)*XB(7)+XB(4)*XB(8) 

F(8)=XB(5)*XB(1)+XB(6)*XB(2)-RB 

RETURN 

END 


Time  Minimization  Vernon 


SUBROUTINE  rCNI(N,T,X,XP) 

DOUBLE  PRECISION  T,X(N),XP(N),MUlRO,B,D1TH,M0,MD,U 
INTEGER  N 

COMMON  /A/  MU  /B/  TH,MO,MD 
RO=DSQRT(X(l)**2+X(2)**2) 

IF  (RO  .EQ.  0)  RO-l.O 
B=TH/(M0-  MD*T) 

U=DATAN2(X(8),X(7)) 

D=X(7)*X(1)  +  X(8)*X(2) 

XP(1)*X(3)*X(9) 

XP(2)=X(4)*X(9) 

XP(3)=(-MU*X(l)/RO**3  +  B*DCOS(U))*X(9) 
XP(4)-(-MU*X(2)/RO**3  +  B*DSIN(U))*X(9) 
XP(5)=»MU*(X(7)/RO**3  -  3.0d+00*X(l)*D/RO**5)*X(9) 
XP(0)«MU*(X(8)/RO**3  -  3.0d+00*X(2)*D/RO**5)*X(9) 
XP(7)=»-X(5)*X(9) 

XP(8)— X(8)*X(9) 

XP(9)«0.0 

RETURN 

END 
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fVMODTlNl  rCNJ(N,T,X,PD) 
DOUBLE  PRECISION  T,X(N),I 


DOUBLE  PRECISION  T,X(N),PD(N 
•  MD.E,  F.G.H.U 
INTEGER  N 

T1MMON  /A/  MU  /B/  TH.MO.MD 
RO-DSQRT(X(l)**2+X(2)**2) 

(F  (RO  EQ  0)  RO— 1.0 
B-TH/jMOMD-'n 
’  -DATAN2(X(8),X(7)) 

F-X(  1)/X(7) 

-(1  0i+00  ♦  F**2) 

D-X(7)‘X(1)  +  X(8)*X(2) 
t>X<8)*X(l)  f  X(7)*X(2) 
;-1<M*00*X(7)*X(l)  +  X(8)*X(2) 
H-1<M-00*X(2)*X(8)  +  X(1)*X(7) 

PHI  l)-00 
PCX  I  21-0  0 
PtXl  1)-X(9) 

PCX  I  4)-0  0 
PCX1  51-0  0 
PH 1  •)— 0  0 
PCX1  T)-00 
PCX  I  81-0  0 
PCX  I  91-X(3) 


T1X(N),PD(N,N),MU,RO,B,C,D,TH,M0, 


pcxa 
PH2 
PCX2 
PCX2 
PCX  3 
PCX  2 
PH  2 
PH3 
PH3 


l)-0  0 
21-0  0 
D-0  0 
4>-X{9) 
5)-0  0 
•1-0  0 
T)-0  0 
•1-0  0 
9>-X<4> 


PTxi.ij-MU*(  10d+00/RO**3  +  3.0d+00*X(l)**2 
•  RO**5)*X(9) 

PH 2  2>-10d*00*MU*X(l)*X(2)/RO**5*X(9) 

t’CXl  11-0  0 

PHI  4)-0  0 

PHI  51-0  0 

PHI  81—0  0 

phi  n-B-DsiNiui-xw/ic^xm^a^xw 
PHI  •)-  B*DSIN(U)/(C*X(7))*X(9) 
phi  91-  MU*X(l)/RO**3  +  B*DCOS(U) 


PD(4,1)«PD<3,3) 

PD<4,2)-MU*(-1.0d+00/RO**3  4-  3.0d+00»X(2)**2 
/RO**5)*X(9) 

PD(4,3)«0.0 

PD(4,4)-0.0 

PD<4,5)-0.0 

PD(4,fl)-0.0 

PD(4,7)— B*DCOS(U)*X(8)/(C*X(7)**2)*X(9) 
PD(4,8)»B*DCOS(U)/(C*X(7))*X(9) 

PD(4,9)— MU*X(2)/RO**3  4-  B*DSIN(U) 

PD(5,l)»MU*(-3.0d+00*G/RO**5+15.0d4-00*X(l)**2*D 
I  /RO**7)*X(9) 

PD(5,2)-MU*(-3.0d+00*E/RO**5+15.0d+00*X(l)*X(2)*D 
I  /RO**7)*X(9) 

PD(5,3)-0.0 

PD(5,4)=0.0 

PD(5,5)«0.0 

PD(5,6)*0.0 

PD(5,7)-MU*(1.0d+00/RO**3  -  3.0d+00*X(l)**2 
\  /RO**5)*X(9) 

PD(5,8)--3.0d+00*MU*X(l)*X(2)/RO»*5  •  X(9) 
PD(5,9)=MU*(X(7)/RO**3  -  3.0d+00*X(l)*D/RO**5) 

PD(«,1)*PD(5,2) 

PD(0,2)-MU*(-3X)d+OO*H/RO**5+15.Od+OO*X(2)**2*D 
I  /RO**7)*X(9) 

PD(8,3)-0.0 

PD(G,4)=»0.0 

PD(8,5)=0.0 

PD(8,«)«0.0 

PD(8,7)=PD(5,8) 

PD(8,8)»MU*(1.0d+00/RO**3  -  3.0d+00*X(2)**2 
\  /RO**5)*X(9) 

PD(8,9)»MU*(X(8)/RO**3  -  3.0d+00*X(2)*D/RO**5) 

PD(7,1)*0.0 
PD(7,2)«0.0 
PD(7,3)»0.0 
PD(7,4)-0.0 
PD(7,5)«- X(9) 

PD(7,«)-0.0 
PD(7,7)-0.0 
PD(7,8)-0.0 
PD(7,9)— X(5) 


PDd.U-00 

PD<8,2)-00 

PD<8.3)-00 

PD<8,4)-00 

PD<8,5)-0.0 

PD(8,e)--X(») 

PD(8,7)-0.0 

PD(8,8)-00 

PD<8,t)— X<6) 

PD(9,l)-0.0 

PD<9,2)-00 

PD<9,3)-00 

PD<9,4)-0.0 

PD<9,5)-0.0 

PD(9,6)-00 

PEX9,7)-0.0 

PD<9,8)-00 

PD(9,9)-0.0 

RETURN 

END 


SUBROUTINE  FCNB(N,XA,XB,F) 

DOUBLE  PRECISION  XAINJ.XBINLFINJ.KS.KM.RB.U.B.XPa, 
8  XP4,TH,M0,MD,MU 
INTEGER  N 

COMMON  /A/  MU  /CONV/  KM.KS  /B/  TH,MO,MD 
RB— DSQRT(XB(1)**2+XB(2)**2) 

U-DATAN2(XB(8),XB(7)) 

B— TH/(M0-MD*XB(9)) 

XP3— XB(9)*XB(7)*(-MU*XB(l)/RB**3+B*DCOS(U)) 
XP4-XB<9)*XB<8)*(-MU*XB(2)/RB**3+B*DSIN(U)) 
F(I)-XA(1)  -  KM 
F(2)-XA(2) 

F(3)-XA(3) 

F(4)-XA(4)  -  KS 

F(5)-XB  1)*XB(3)+XB(2)*XB(4) 

F(8)-XB(7)*XB(3)+XB(8)*XB(4) 

F(7)-RB  -  6.701 18388d+00 

F(8)-XB(5)*XB(1)+XB(6)*XB(2)  -  RB 

F(9)-XB(9)+XB(9)*(XB(5)*XB(3)+XB(6)*XB(4))+XP3+XP4 

RETURN 

END 


1 
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SUBROUTINE  K£PLER(T,MU,RO,VO,R,V) 

DOUBLE  PRECISION  T,R0(4),V0(4),R(4),V(4),E,MU, 
t  SMU,A,XO,XN,Z,SZ1C,S1TN1DT,DOTITOL1F,G)FD,GD 
if  (T  eq.  0)  then 
R(l)-RCHl) 

R(2)«RO(3) 

R(3)-R0(3) 

V(1)-V0(1) 

V(2)»V0(3) 

V(3)«V0(3) 

R(4)«DSQRT(R(1)**2+R(2)**2+R(3)**2) 

RETURN 

endif 

DOT-RO(1)*VO(1)+RO(2)*VO(2)+RO(3)*VO(3) 
E-VO(4)**2/2.0IHOO  -  MU/R0(4) 

A— -  1.0D+00/(2.0D+00*E) 

SMU-DSQRT(MU) 

TOL-l.OD-lO 
XO-SMU  •  T/A 
Z»X0**2/A 
SZ-DSQRT(Z) 

C— (l.d+00  -  DCOS(SZ))/Z 
S-(SZ  -  DSIN(SZ))/SZ**3 

TN-DOT/SMU*XO**2*C  +  (l.d+0O-R0(4)/A)*XO**3*S 
I  +  R0(4)*X0 

DT-X0**2*C  +  DOT/SMU*X0*(l.d+00-Z*S)  +  R0(4) 

I  •  (l.d+00-Z*C) 

XN-XO  +  (T-TN)/DT 
IF  ((XN-XO)  GT.  TOL)  THEN 
XO-XN 
GO  TO  10 
ENDIF 

F-l.d+00  -  XN**2/R0(4)  •  C 

G«XN”2  •  DOT/MU  •  C  +  R0(4)*XN  •  (l.d+00-Z*S)/SMU 
DO  20  1-1,3 

R(I)-F  •  R0(I)  +  G  •  V0(I) 
R(4)-DSQRT(R(1)**2+R(2)**2+R(3)**2) 

GD— l.d+00  -  XN«2/R(4)  •  C 

FD— SMU/(RO(4)*R(4))  •  XN  •  (Z*S- l.d+00) 

DO  30  1-1,3 

V(I)-FD  •  R0(I)  +  GD  •  V0(I) 

RETURN 

END 
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